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conditions  at  "artificial^  boundaries  for  semi-linear  elliptic  boundary  valueN'-— ^ 
problems  on  semi-infinite  cylindrical  domains.  A  general  theory  developed  by 

the  authors  in  a  previous  work  [11]  is  applied  to  establish  the  existence  of 
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exact  boundary  conditions  and  to  obtain  useful  approximations  to  them.  -They 
are  based  on  the  Laplace  transform  solution  of  the  linearized  problem  at  in- 
f inity.  -We  discuss  the  incorporation  of  these  conditions  in  a  finite  differ¬ 
ence  scheme  and  present  the  results  of  a  numerical  experiments  the  solution 
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of  the  Bratu  problem  in  a  two  dimensional  stepped  channel.  We  also  examine 
certain  problems  concerning  the  existence  of  solutions  on  infinite  domains. 
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SIGNIFICANCE  AND  EXPLANATION 


Many  computational  problems  arising  in  applied  mathematics  are  posed  in 
infinite  pipes  and  channels.  One  of  the  most  important  examples  of  this  is 
the  problem  of  incompressible  fluid  flow  in  such  a  geometry.  As  computations 
are  only  possible  on  finite  domains,  "artificial"  boundaries  must  be  intro¬ 
duced  and  boundary  conditions  must  be  imposed  there.  For  the  fluids  problem 
these  are  referred  to  as  inflow  and  outflow  boundary  conditions. 

In  a  previous  work  the  authors  developed  a  general  theory  of  boundary 
conditions  at  artificial  boundaries.  In  this  work  we  show  how  to  apply  that 
theory  to  the  numerical  solution  of  semi-linear  elliptic  problems.  Such  prob¬ 
lems  are  well-suited  for  numerical  experimentation  for  a  variety  of  reasons: 
first,  the  abstract  theory  can  be  directly  applied  to  them;  second,  the 
derivation  of  boundary  conditions  for  these  problems  is  formally  applicable  to 
the  equations  of  incompressible  flow;  third,  the  problems  are  physically  and 
mathematically  interesting  in  their  own  right. 

We  illustrate  the  large  reduction  in  the  error  brought  about  by  use  of 
the  asymptotic  boundary  conditions  through  presentation  of  the  results  of 
computations  on  the  Bratu  problem.  This  deals  with  the  existence  of  stable 
solutions  to  the  equation  V^u  *  ~Xeu,  which  is  used  as  a  model  of  thermal 
ignition.  For  X  sufficiently  large  and  positive,  solutions  do  not  exist 
and,  hence,  ignition  is  said  to  have  taken  place.  The  problem  is  to  determine 
the  minimum  value  of  X  for  which  this  occurs.  We  present  some  results 
concerning  the  change  in  the  critical  value  of  X  due  to  finite  perturbations 
of  infinite  domains. 


The  responsibility  for  the  wording  and  views  expressed  in  this  descriptive 
summary  lies  with  MFC,  and  not  with  the  authors  of  this  report. 


NUMERICAL  SOLUTION  OF  SEMI-LINEAR  ELLIPTIC 
PROBLEMS  ON  UNBOUNDED  DOMAINS 


Thomas  M.  Hagstrom*  and  H.  B.  Keller** 


1.  Introduction 


We  consider  boundary  value  problems  of  the  form: 


(1.1) 


a)  Lu( x,y )  =  f(u,£),  (x»x)  e  [0,»)xQ,  ft  C  Rn“  ; 

b)  aft(X}  7$  +  bn(x)u  =  X  e 

c)  ao(*}  17  +  V*)u  "  C0(X) •  x  ■  °» 

. .  lim  ,  ,  .  .  lim  3u  ,  .  A 

d>  u(x'X}  =  U»(X>'  .  7x  l*'X]  =  0  * 


Here,  L  is  a  linear  uniformly  elliptic  second-order  operator  which  is 
independent  of  x: 

*’  L  5  72 +  L.  •  ♦  V 


(1*2>  b>  L1  "  Jiani(X)a?7  +  an(^)? 


n-1  n-1  .2  n-1  a 

l2  1  l  i.  +  l  ai(x)^r  +  a(x> 

i= 1  j=1  J  yi  yj  i=l  1  oyi 
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The  expression  denotes  the  conormal  derivative  associated  with  L  and 

ft.  The  function  u^fy)  must  satisfy  the  limiting  cross-sectional  problem 
obtained  by  formally  letting  x  00  and  setting  u^  =  u^^  20  in  (1.1a,b)  to 
get : 

a)  l2u®  =  X  e  a' 

d.3)  aUoe 

b)  aft(^)  a TT  +  bft(x)u«  =  cn(X>'  x  e  aft  . 

Equations  such  as  (1.1a)  arise,  for  example,  in  equilibrium  problems  in  non¬ 
linear  heat  generation.  (See,  e.g.,  Aris  [3].)  Furthermore,  the  method  we 
describe  can  be  generalized  to  systems  and  to  higher  order  equations  [11].  As 
such  it  has  potential  for  application  to  steady  state  problems  of  fluid  flow 
in  pipes  and  channels. 

Our  procedure  is  to  introduce  an  artificial  boundary  at  some  point 
T  >  0,  impose  boundary  conditions  there  and  solve  the  resulting  finite 
problem  on  [0,t]  x  ft  by  a  standard  numerical  technique.  A  theory  of  boun¬ 
dary  conditions  at  an  artificial  boundary  has  been  developed  by  the  authors  in 
[11]  and  we  use  it  here.  Indeed,  one  of  the  purposes  of  this  work  is  to 
illustrate  the  power  of  the  general  theory. 

Other  authors  have  discussed  the  problem  of  deriving  boundary  conditions 
at  an  artificial  boundary  for  linear  elliptic  problems.  Gustafsson  and  Kreiss 
[10]  point  out  the  possibility  of  deriving  exact  conditions  by  use  of  a 
Laplace  transformation  in  x.  Fix  and  Marin  [7]  and  Goldstein  [9]  use  a 
related  approach  to  solve  problems  in  underwater  acoustics  and  wave  propa¬ 
gation  in  cylindrical  waveguides.  The  first  approximation  to  the  boundary 
condition  given  by  our  method  coincides  with  that  obtained  by  applying  the 
Laplace  transform  method  to  the  problem  linearized  about  u^fy). 


In  section  2  of  this  work  we  state  the  basic  results  of  [11]  and  show  how 
they  can  be  applied  to  (1.1).  In  particular,  we  present  detailed  asymptotic 
expansions  of  the  exact  boundary  conditions  and  state  sufficient  conditions 
for  their  validity.  Extensions  to  other  problems  are  also  discuised.  In 
section  3  the  inclusion  of  the  expansions  in  a  discrete  approximation  is 
described  and  their  efficient  use  is  considered. 

A  special  case  of  (1.1a),  the  Bratu  problem  in  a  two-dimensional  channel, 
is  introduced  in  section  4.  We  discuss  its  physical  interpretation  and  quote 
various  existence  results  for  finite  domains.  The  results  of  some  numerical 
experiments  are  presented  assessing  the  effects  of  varying  the  location  of  the 
artificial  boundary  and  the  number  of  terms  in  the  asymptotic  expansion.  Some 
questions  concerning  the  existence  of  solutions  in  unbounded  domains  are  also 
examined. 


2.  Construction  of  the  Boundary  Conditions 

In  order  to  conform  to  the  notation  of  [11],  we  write  the  problem  as  a 

first  order  ordinary  differential  equation  in  a  Banach  space.  Letting 

v  =  u  -  u„  we  introduce: 

00 


|!f(x'xA  =  /wi  \ 

\v(x,y)  /  \w2  / 


We  seek  a  solution,  v,  which  for  each  x  is  an  element  of  the  space 
w°(ft)  -  that  subspace  of  the  Sobolev  space  W2(ft)  consisting  of  functions 
which  satisfy  the  homogeneous  version  of  (1.1b),  i.e.  with  =  0.  Here  for 
each  x  we  require: 


w 

a)  w( x )  e  8  =  {C^1 ) :  w  e  W  (0),  w2  e  W2(0),  (2.2b)  holds.} 


(2.2) 


n-1  n-1  3w2  n-1 

b)  a0  V  V  n.a.  .  -s —  +  an  V  n.a  ,w.  +  b0w_  =  0,  y  e  3ft 
ft> „  ,L .  i  it  3y .  .  l  ni  1  ft  2  * 

i=1  3=1  3  i=1 


Here  nA  is  the  ith  component  of  the  unit  normal  to  3ft.  Choosing  t  >  0  as 
the  location  of  the  artificial  boundary,  we  rewrite  the  equation  in  the  tail 
in  the  abstract  form: 


(2.3) 


dw 

a)  —  =  Aw  +  P ( w ) ,  x  >  T ; 


b)  ^  w(x)  =  0; 

X-H» 


where 


a)  A 


=  rL1W1  -  L2W2  +  fu(U-)W2\  . 

w2  \  wi  / 


(2.4) 


/f(u_  +  w_)  -  f  (u  )w„  -  f(u  )\ 


In  [11]  an  explicit  approximation  to  the  solution  of  (2.3)  is  constructed 


in  terms  of  the  eigenfunctions  of  A.  That  is,  we  consider  pairs 
(Ajj,  wa)  e  <t  x  8  such  that: 

(2.5)  Aw^  =  A^w i  =  1,2,...;  11  U  *  0  . 

We  assume  that: 

a)  The  eigenvalues,  A^,  are  distinct  and  bounded  away  from  the 
imaginary  axis. 

(2.6) 

b)  The  eigenfunctions,  w^,  form  a  Riesz  basis  of  8.  (See, 
e.g.,  Gohberg  and  Krein  [8]  for  a  discussion  of  non-orthonormal 
bases . ) 


Note  that  (2.6a)  ensures  that  the  linearized  problem  has  an  exponential 
dichotomy.  From  (2.4a)  we  see  that  (2.5)  is  equivalent  to  the  eigenvalue 
problem  obtained  by  application  of  a  Laplace  transformation  in  x  to  the 
linearization  of  (1.1)  about  u  : 

CO 


V*  '  VU»)Y*  +  VlY£  +  V*  =  °'  *  e 


n- 1  n- 1  3  ] 


(2.7)  b 


1  ni<(  l  aij  af5  +  VniY*t  +  Vi  •  »•  X  « 

1=1  3=1  3 


c)  W. 


Vt 


Condition  (2.6b)  can  be  particularly  difficult  to  check  for 
nonself ad joint  problems.  (See,  e.g.,  Berezanskii  [6]  for  a  discussion  of  the 
selfadjoint  case.)  Agmon  [1]  treats  the  case  =  0  while  Agmon  and 

Nirenberg  [2]  give  sufficient  conditions  for  completeness  in  the  class  of 


initial  data  leading  to  absolutely  integrable  solutions  in  the  tail.  Note 


that  the  set  of  all  decaying  solutions  of  the  linearized  problem  is  generated 


by  those  w^  whose  eigenvalues  have  negative  real  part.  We  denote  this 
subspace  of  8  by  A  : 

(2.8)  A  =  spantw^:  ReA^  <  0}. 

We  now  state  our  fundamental  result  concerning  the  existence  of  solutions 
of  (2.3): 

Theorem  2.9  Let  6  >  0  satisfy 


(2.9)  6  <  - - - - -  , 

*  S  +  K3  >(-*-) 

+ 

where 


a ) 


X+  = 


inf 

Re  A^>0 


RefA^) , 


A  = 


inf 

ReA£<0 


( -Re ( A . ) ) ; 


b’  $  2|?|<2TS  ’ 


n-1 

+  v,X)(  ; 

c) 

5  i, 

^SUP 

|v | <2y6 

|fuuy.(u- 

n-1 

d) 

K3  5  l 

sup  1 

(».<);>) 

,~S,UP  ,  |fuuu(u~  +  XM 

i=1 

x  e  fi  1 

yi 

1  |v|OY«  1  1 

e)  Y  is  a  constant  appearing  in  Sobolev's  inequality  (see  Agmon  [1]), 


M  <  ^Bu"w2(S2)  + 


"uV(flJ  ? 


f)  a. ,  a  are  constants  such  that  if  w  =  )  c.w, 
12  %  1  1 


then : 


(Such  a  6  clearly  exists  if  the  necessary  derivatives  of  f  are  continuous 
near  u,,,.)  Then,  for  any  E,  e  A  which  also  satisfies: 

(2.11)  “U8  <  IS  ' 

there  exists  a  solution,  w(x),  of  (2.3)  such  that,  if  is  the  projection 

operator  into  A  whose  nullspace  is  the  span  of  the  w^  whose  eigenvalues 
have  positive  real  part,  then: 

(2.12)  Q^wd)  =  C. 

Furthermore,  w(x)  is  unique  among  small  solutions  satisfying  (2.12).  / 

We  note  that  Theorem  2.9  is  the  specialization  of  some  results  of  [11]  to 
the  present  case.  The  proof,  as  given  there,  follows  this  path:  [i)  a 

solution,  wq(x),  of  the  linearized  problem  satisfying  wQ(x)  =  5  is  written 

down,  (ii)  w(x)  is  represented  as  wQ(x)  +  w(x),  (iii)  an  integral  equation 
for  w  is  derived  from  an  integral  representation  (in  this  case  a  Green's 
function  representation)  of  solutions  of  the  linearized  inhomogeneous  equations, 
(iv)  the  existence  of  w  is  established  by  a  contraction  mapping  argument. 

In  order  to  carry  out  step  (iv),  it  is  necessary  to  prove  certain  estimates. 

This  leads  to  the  inequality  (2.9).  As  their  derivation  is  straightforward, 
we  postpone  its  presentation  to  the  appendix. 

The  solution  we  have  constructed,  w(xj£),  can  now  be  used  to  write  down 

an  exact  boundary  condition  at  x  =  x.  Written  in  abstract  form  it  is 

( see  [11]  ) : 

(2.13)  (I-Pb)w(t)  =  -/"  S(T,p)  (I-QJ  r(w(P,Qoow(t)  ))dp, 

Here,  the  combination  Sfs.tMl-Q,,,)  is  the  solution  operator  of  the 
linearized  problem  restricted  to  the  span  of  w^  whose  eigenvalues  have 
positive  real  part.  Its  existence  is  guaranted  only  for  s  <  t.  Equation 


I 


1 


I 

t 

r 

r 


* 


(2.13)  is  an  exact  condition  in  that  there  exists  a  small  solution,  w,  of 
(2.3)  if  and  only  if  w(t)  satisfies  (2.13). 

A  useful  approximation  to  (2.13)  is  obtained  in  terms  of  the  eigen¬ 
functions,  .  The  iterative  process  implicit  in  the  contraction  argument  is 
truncated  to  give  some  approximation,  v/n^(x;£),  to  w(x;£).  Furthermore, 
the  nonlinearity  R  is  replaced  by  a  finite  Taylor  series.  Then  the 
integrals  involved  only  contain  exponential  functions  of  the  integration 
variable  and  can  be  evaluated  explicitly.  The  result  is  an  expression 
relating  the  expansion  coefficients.  (See  [11],  eq.  (6.1  .  e  rewrite  the 

quadratic  approximation  below  for  an  expansion  in  terms  ..  Y#(^).  This,  in 

turn,  makes  use  of  solutions  of  the  adjoint  problem: 


*  _ 


a)  L2  Yi  "  fu(u~)Y£  +  X£  L1  Yi  +  X7yi  =  X  e  Q; 


(2.14) 


3Y„ 


n-1  n-1  3a. 


b)  an  +  (ba  +  an  E  W1  "  ai^  5i  =  °'  * 

i=1  3=1  3 


e  a n. 


Here,  Lj  and  are  the  formal  adjoints  of  L2  and  Li*  We  choose  our 

normalization  so  that: 

n- 1 


a)  /  dy[(X  +  A  +  a  )Y  -  | — (a  .  Y.  )  ]  Y  =5  ; 

*  **  £  m  n  £  3y  mi  m  £m 

ii  1=  1  1 


(2.15) 


b)  II Y  II,.  ,0,  =  1. 
m  W2  (ii) 


(VX(T^>\ 
W(T,y)  } 


Given  the  pair  \  "t  ~  x  j  we  define  expansion  coefficients  in  the  following 
way  : 

c5  =  /  dX 

x  a 


(2.16) 


>-*  n~1  3 


+  i +  VxOVxJ  -  i  j^r(ani(x>Y1(x))]v(‘t'X>l  • 

i=1  y  i 


-8- 


1 


We  further  define  matrix  elements  of  the  quadratic  approximation  to  R: 


(2.17) 


=  4-  /  dvY  (y)  f  (u  ,y)Y  (y)Y  (y) 

mn  2  ^  l  uu  “  ^  m  ^  n  ^ 


a 


We  then  have  as  an  0  ( II  v(  T  ,^ )  11  ^ )  approximation  to  (2.13): 


1 


t 


L 

i 


(2.18) 


( £ )  CmCn 

ci  =  i  l  a  mn  x  +x  -x.  '  v£  such  that  ReXe  >  0l 

m  n  m  n  l 

ReX  <0  ReX  <0 
m  n 


Equations  (2.16)  and  (2.18)  take  on  a  more  familiar  form  when  the 
eigenvalue  porblem,  (2.7),  is  selfadjoint.  In  particular  we  consider  the  case 
when  L  is  the  laplacian  and  the  boundary  conditions  are  of  Dirichlet  type. 

We  then  have: 

n- 1  32Y 

al  i  —  -  +  V*  -  "■  xe  “> 


i=  1  ■'i 


(2.19) 


t>)  y  (X)  =0,  x  e  an. 


Define  <f>£  =  y£Y£  to  be  normalized  in  L„  ft) 


(2.20  ) 


I  dX  =  K 

ft  x 


Then  we  have  by  (2.15a): 


V*1  ■  *i  ■ 


(2.21) 

Note  that  each  eigenfunction  Y£  corresponds  to  two  distinct  eigenvalues, 

v 

±X  .  Consider  the  expansion  of  the  pair  (  )  in  terms  of  the  two  bases: 


(2.22) 


/v  (x,y)\  ,  ^  Y 

U,x, )  =  itKc«)(  i  Vv*)( 


rh£(x), 

Y£  ''  '  £  ^£<x)J  ** 


The  expansion  coefficients  are  related  by: 


-9- 


4 


I 

I 


1 

» , 
L  ' 
►  ' 
k 
. 

I 

l" 


(2.23) 


Define  matrix  elements  in  the  new  basis  by: 


(2.24)  5^  =  4  /  dy  f  (u  ,y)  <j>.(v)  4>  (y)  <)>  (y). 

mn  2  '  ^  uu  00  *■  Z  ^  m  *•  n  'L 


They  are  related  to  the  old  matrix  elements  through: 

y  t 


(2.25) 


*  s(t) 


,<*>  _ _ 

mn  2X„Y  Y  ^mn 
x.  m  n 


Using  (2.23)  and  (2.25)  the  boundary  condition  (2.18)  becomes: 


(2.26) 


S  +  V*  =  -I  l 


-(£) 

a 

mn 


X  +X  +X .  4^“  m  X 
m  n  m  n  l  m 


-  b1  h* 

U*  -  T1)  (h  -  T“ )  , 


which  can  be  replaced,  to  the  same  quadratic  order  of  approximation,  by: 


(2.27)  h^  =  -X£h£  l  l  X  +X  +1  hmhn 

m  n  m  n  l 


Condition  (2.27)  is  used  in  the  calculations  presented  in  section  4. 

In  closing  this  section  we  note  that  the  analysis  given  above  can  easily 
be  extended  to  a  variety  of  different  situations.  The  assumption  concerning 
the  lack  of  dependence  of  the  domain  and  coefficients  on  x  is  necessary  only 
in  the  tail,  i.e.  for  x  sufficiently  large.  In  fact,  the  asymptotic 
expansions  given  in  [11]  allow  for  a  dependence  of  the  coefficients  on  x  so 
lonq  as  they  approach  their  limiting  values  sufficiently  fast.  Also,  the 
restriction  to  scalar  equations  in  not  needed. 
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3.  Discrete  Approximations 

We  show  how  to  implement  the  boundary  conditions  discussed  above  in  a 
numerical  computation.  For  simplicity  we  assume  a  one  dimensional  cross- 
section,  ft,  with  Dirichlet  boundary  conditions.  We  also  consider  centered 
finite  difference  approximations  to  the  derivatives  with  a  uniform  mesh 
width.  Let  p  denote  the  number  of  gridpoints  in  a  cross-section,  h  be  the 
mesh  width  and  i,  1  <  i  <  p,  parametrize  the  points  in  a  cross-section. 
Then,  if  j  parametrizes  the  x  coordinate,  the  index,  k,  of  a  mesh  point 
is  given  by  k  =  p( j— 1 )  +  i.  At  an  interior  point  a  finite  difference 
approximation  to  (1.1a)  is  given  by: 


1  r  s  a2(yi) 

— (u(k+p)  -  2u(k)  +  u{ k-p) )  + 


2h 


ai . (y. ) 

+  - - - (u(k+1 )  -  2u( k )  +  u(k-1)) 


(u(k+p)  -  u(k-p) ) 


(3.1) 


a21 (yi - 


+  h2 


a1(yi) 


(u(k+p+1 )  +  u(k-p-l)  -  u(k-p+1 )  -  u( k-p+ 1 ) ) 


2h 


(u(k+1)  -  u(k-1 ) )  +  aty^Jufk)  -  f(u(k),y  )  =  0 


In  order  to  implement  (2.18),  it  is  necessary  to  solve  the  eigenvalue 
problems  (2.7)  and  (2.14).  We  approximate  these  on  the  same  cross-sectional 
mesh  as  the  full  equation.  We  must  therefore  assume  that  the  eigenvalues  and 
eigenfunctions  which  are  approximated  well  on  this  mesh  are  enough  to  resolve 
the  solution.  (See  Kreiss  [13]  for  a  discussion  of  the  approximate  solution 
of  nonselfad joint  eigenvalue  problems.)  We  denote  our  approximate  solutions 

A  A  A  ^  A 

to  (2.7)  and  (2.14)  by  the  pairs  ( X ^ ,  Y^fy^)),  (\*  ^£*yi^'  Approximate 

expansion  coefficients  are  given  by: 
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i 


(3.2) 


c*  =  F  ^  {^(yi)vx(T'yi)  +  +  a2(yi^  Vyi> 

i=1 

-  k  ^a2i(yi+i}  VW  -  a2i(yi-i}  Vyi-i,)lv(T'yi)^ 


^  (  il )  (  £  ) 

and,  if  approximations,  or  ,  to  the  matrix  elements  ot^  are  calculated 
in  a  similar  fashion,  (2.18)  can  be  replaced  by: 


(3.3) 


Ct  =  1 


V  >  aU)  ~\n  -  ,  Vi  such  that  ReA  >  0. 

®n  w  0 

m  n  A  + A  -A. 

,  ,  m  n  x. 

ReA  <0  ReX  <0 
m  n 


We  note  that  v(T,y^)  =  u(T,y^)  -  u^fy^)  and  that  vx(T,y^)  must  be  replaced 
by  a  difference  approximation.  A  simple  approach  is  to  take  T  half  way 
between  vertical  grid  lines  and  to  calculate  v(x,y^)  by  averaging.  In 
general,  there  will  be  2p  eigenvalues,  A^,  of  the  discrete  problem. 
Equation  (3.3)  represents  as  many  equations  as  there  are  eigenvalues  with 
positive  real  part.  If  this  is  not  equal  to  p,  then  we  cannot  expect  the 
discrete  equations  to  be  well-posed.  (One  reason  for  this  might  be  that  the 
original  problem  is  ill-posed. ) 

As  expected,  all  of  this  is  simplified  in  the  selfadjoint  case.  Then, 
(3.3)  can  be  replaced  by  a  discrete  analogue  of  (2.27)  where 


h'i  =  ?  Vyi,Vx(T'yi)j 

i=  1 


(3.4) 


h*  =  ?  X  n(yi)v(T'yi)- 

1=1 


Since  the  p-vectors  ^(y^)  are  orthogonal,  our  approximation  to  (2.27)  can 
be  rewritten  in  the  convenient  form: 
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(3.5) 


v 

X 


(T,yi) 


-  f  1  My,  )M4(y-i)v(T,yi) 

4-1  j-1  1  *  *  3  3 


l  l  f 

4=1  m= 1  n=1 


a  d  p 

My,)-* — T”  a  (  l  4>  (y.  )v(x,y  ))(  l  $  (y  )v(T,y  )j  . 

1  1  X  +A  +X„  j-1  m  3  3  j-1  n  3  3 

m  n  4  J  J 


We  now  consider  the  solution  of  the  system  of  nonlinear  equations  by 
Newton's  method.  From  (3.1)  we  see  that  the  interior  of  the  Jacobian  will  be 
banded  with  band  width  at  most  p  +  1.  Assuming  we  use  (3.3)  to  relate  two 
vertical  grid  lines,  the  last  p  rows  will  have  nonzero  elements  in  the 
last  2p  columns,  increasing  the  band  width  by  as  much  as  p.  If,  however, 
(3.5)  can  be  used,  it  is  possible  to  write  the  condition  in  s.  ch  a  way  that 
the  band  width  is  not  increased.  Therefore,  in  the  latter  case,  no  extra  work 
is  needed  to  solve  the  system  by  banded  Gaussian  elimination.  In  the  former, 
a  bordering  technique  is  necessary  to  avoid  significant  additional  calcula¬ 
tion.  The  effect  of  the  nonlocal  boundary  conditions  on  the  performance  of 
iterative  techniques  has  not  been  examined  in  general.  Bayliss,  Goldstein  and 
Turkel  [5]  have  found  that  the  use  of  the  linear  version  of  (3.5)  has  essen¬ 
tially  no  effect  on  the  convergence  of  their  preconditioned  conjugate  gradient 
algorithm  for  the  Helmholtz  equation  in  three  dimensions. 

A  significant  number  of  new  calculations  are,  however,  required  to 
evaluate  the  Jacobian.  From  (3.5)  we  see  that,  for  the  quadratic 
approximation,  this  involves  an  evaluation  of 


5  J,  *ij»  v"'y«> 


(3.6) 


~(m) 

=  l  1  I  My,)  ^ — *-  My»>  My.) 


m=1  n=1  K= 1 


X  +X  +  X 
m  n  K 


Since  A.  ..  can  be  evaluated  beforehand,  this  entails  O(p^)  new 
1 3  * 

multiplications  at  each  iteration.  The  direct  solve  itself  requires  0(p3q) 
operators  where  q  is  the  number  of  gridpoints  in  the  x  direction.  There 
are  two  ways  to  reduce  the  error  due  to  the  introduction  of  the  artificial 
boundary:  first,  to  take  more  terms  in  the  expansion  of  the  boundary 
condition;  second,  to  move  the  boundary  further  out.  The  latter  requires 
0(p3Aq)  additional  multiplications.  The  increased  cost  of  evaluating  the 
Jacobian  in  thie  former  approach  is  0(ps+^  where  s  is  the  number  of  terms 
taken  in  the  expansion.  This  suggests  that  the  quadratic  approximation  is  an 
efficient  choice.  These  considerations  might  change,  of  course,  if  a 
different  solver  was  used  or  if  the  dimension  of  the  cross-section  were 
higher. 


4.  The  Bratu  Problem 


We  consider  the  following  special  case  of  (1.1a): 


(4.1  ) 


32u  t  32u  _  _Xeu 
3x2  3y2 


Equation  (4.1),  which  is  also  associated  with  the  names  of  Gelfand  and  Frank- 
Kamenetskii ,  arises  in  the  theory  of  thermal  ignition  of  gases.  (See,  e.g., 
Aris  (3].)  The  problem  of  existence  of  positive  solutions  has  been  considered 
by  various  authors.  We  state  below  two  theorems  from  the  literature  for  the 
Dirichlet  problem  on  a  finite  domain  -  that  is  (4.1)  holds  on  some  finite 
domain,  D,  and: 

(4.2)  u  =  0,  (x,y)  e  3d. 

Theorem  4.3  (Keller  and  Cohen  [12])  Let  X  >  0  be  such  that  (4.1,  4.2)  has 

a  positive  solution.  (We  say  that  X  is  in  the  spectrum.)  Then,  if 
*  * 

0  <  X  <  X,  X  is  in  the  spectrum.  Furthermore,  for  all  X  in  the 
spectrum,  there  exists  a  minimal  positive  solution,  uQ(x,yf  X),  such  that, 
if  U(x,yj  X)  is  any  other  positive  solution,  then 

(4.3)  Uq (x,y;  X)  <  U(x,y;  X)  V(x,y)  B  D. 


The  minimal  positive  solution  is  stable  in  that  the  linearized  eigenvalue 


problem 

a)  V2<Ji  +  Xeu<j>  =  a$  in  D; 

(4.4) 

b)  $  =  0  on  3D; 


has  only  negative  eigenvalues. 

Theorem  4.5  (Bandle  (4))  Let  D'  £  D.  Then,  if  X  is  in  the  spectrum  for 
D,  it  is  in  the  spectrum  for  D'.  Furthermore,  the  minimal  positive 
solutions  satisfy: 

(4.5)  un(x,y;  X,D')  <  u»(x,y;  X,D)  V(x,y)  e  D' . 


That  is,  the  minimal  equilibrium  temperature  increases  with  the  size  of  the 
doma i n . 

We  consider  the  problem  (4.1)  on  domains  with  infinite  boundaries  of  the 
following  types: 


(4.6) 


ye  (0,1),  |x  -  xn I  >  a;  ye  [d, 1 ) ,  |x  -  xn|  <  a; 


d  >  0 


d  <  0 


where  d  can  be  positive  or  negative.  On  the  top  and  bottom  of  the  channel 
and  on  the  step  at  x  =  x„  +  a  we  require  u=0.  As  [  x  [  -*•  «»  we  require: 


(4.7) 


u(x,y)  =  u  (y)  . 


Here,  uM(y)  is  a  solution  of  the  limiting  cross-sectional  problem: 


(4.8) 


a )  u"  =  -Ae  ,  y  e  ( 0 , 1 ) ; 


b)  u  (0)  =  u_(1 )  =  0. 


We  are  now  assuming  that  A  is  in  the  spectrum  of  the  cross-sectional 
problem  at  infinity.  Furthermore,  in  order  to  ensure  that  condition  (2.6a)  is 


met,  we  must  take  u^,  to  be  the  minimal  positive  solution.  It  is  possible  to 
solve  (4.8)  analytically.  The  spectrum  is  found  to  range  between  0  and  Ac 
where 


(4.9) 


A  =  3.51  ...  . 

c 


.0.- 
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Finally,  we  seek  solutions  which  are  symmetric  with  respect  to  xQ.  This 
allows  us  to  consider  solutions  on  the  semi-infinite  domain  created  by  the 
restriction  of  the  original  domain  to  x  >  Xg.  At  Xg  we  impose  the  boundary 
condition: 

(4.1°)  |£<Vy)  =  °- 

For  the  numerical  solution  of  the  problem  an  artificial  boundary  is 
introduced  at  the  point  x  =  T.  Three  different  boundary  conditions  are 
imposed  there : 

(i)  the  quadratic  condition,  (2.26); 

(ii)  the  linear  approximation,  h^  +  X^h^  =  0; 

(iii)  the  "naive"  zero-order  condition,  =  0. 

Before  discretization,  the  stepped  channel  was  mapped  to  a  straight 
channel  using  the  conformal  mapping,  with  d  >  0: 

a)  s+it=w=F1(z),  z  *  x  +  iy; 

b)  F(w)  =  { £n[  ( e7Tw-1 ) 1//2  +  (enw-  ( 1-a ) ) 1//2  ]  - 

(4.11) 

,,  ,1/2„  r  (  it  w  ,  ,'>1/2  .  ,,  ,1/2,  ttw  ,.1/211  ,  ,,  ,1/2 
(1-a)  '  An  [[e  -(1-a)J  +  (1-a)  (e  -1)  Jj  +  (1-a)  w 

c)  a  =  2d  -  d2  . 

For  d  <  0,  set  d*  =  ■  and  replace  F  by  - — -  pconjg( _wconi<3;  d*)  +  d. 
The  straight  boundaries  at  x  =  xQ,  t  are  mapped  into  slightly  curved 
boundaries.  The  resulting  perturbation  of  the  boundary  conditions  are 
calculated  using  linear  interpolation.  The  equations  are  discretized  and 
solved  (Newton's  method  and  Gaussian  elimination)  as  described  in  section  3. 

For  all  cases  shown  the  uniform  mesh  width  is  h  =  .05.  Due  to  the  curvature 
of  the  boundary,  it  is  necessay  to  implement  the  condition  at  x  =  i  in  such 
a  way  that  the  band  width  is  increased.  To  avoid  additional  computation,  a 


bordering  technique  is  used  to  solve  the  Newton  iterates.  All  calculations 
were  performed  on  VAX  11/780  computers  at  the  California  Institute  of 
Technology  and  the  Mathematics  Research  Center. 

If  the  step  parameter,  d,  is  taken  to  be  zero  a  trivial  solution  exists  - 
namely  u  =  u^fy).  Theorem  4.5  (which  has  not  been  established  for  an 
unbounded  domain)  suggests  that  solutions  exist  if  0  <  d  <  1.  We  always 
found  this  to  be  the  case.  Presented  below  are  results  d  =  .4  and  A=  3.51, 

very  close  to  the  critical  value  for  the  cross-section.  The  minimum 

eigenvalue  of  the  cross-sectional  problem  is  .66825,  so  the  decay  to  u^  is 

relatively  slow.  In  Table  I  we  list  the  maximum  error  as  a  function  of  T 
and  the  order  of  the  asymptotic  boundary  condition.  (The  exact  solution  is 
taken  to  be  a  finite  difference  solution  on  a  large  domain,  x  =  2.262,  using 
the  quadratic  approximation  to  the  boundary  condition . ) 


X 

#  vertical  gridlines 
to  the  right  of  step 

B.C.  approximation 
order 

max  |  h  | 

X=T  lu  “UJ 

„  h  h  „ 

Hu  -u  1 

approx.  “> 

0:  u*  =  0 

.0743 

1.262 

20 

1:  linear 

.076 

.0030 

2:  quadratic 

.0015 

0 

.  1472 

.764 

10 

1 

.177 

.0082 

2 

.0018 

0 

.1936 

.520 

5 

1 

.269 

.0154 

2 

.0051 

0 

.2131 

.381 

2 

1 

.341 

.0213 

2 

.0105 

TABLE  I 


We  see  that  the  quadratic  condition  is  consistently  the  best  and  that  the 
naive  condition  is  consistently  the  worst.  The  success  of  the  quadratic 
condition  is  graphically  displayed  in  Figures  1-4.  Figure  1  shows  the 
solution  using  the  solution  using  the  quadratic  condition  on  a  1 nrge  domain, 
x  =  1.262.  (The  step  is  located  at  x  =  -.057.)  Plotted  are  level  curves 
of  u.  In  Figure  2  the  same  level  curves  are  plotted  for  a  solution  on  a 
small  domain,  t  =  .381.  Figure  3  shows  the  superposition  of  the  two 
solutions.  As  predicted  by  the  linear  error  analysis  (Hagstrom  and  Keller 
[11]),  the  error  decays  as  we  move  into  the  interior.  Figure  4  shows  the 
superposition  of  the  large  domain  solution  and  a  solution  on  the  small  domain 
found  using  u  (x,y)  =  0.  The  level  curves  are  seen  to  be  greatly  distorted. 


FIGURE  3 


FIGURE  4 


A  more  interesting  case,  from  the  point  of  view  of  existence  theory, 
occurs  when  d  <  0.  Then  we  may  expect  that  the  presence  of  the  finite  thick 
region  will  preclude  the  existence  of  steady  solutions,  even  though  the 
infinite  part  is  thin  enough  for  a  one  dimensional  solution  to  exist 
(A  <  A  ).  In  terms  of  the  physical  problem  this  says  that  a  finite  area  of 
excess  thickness  can  cause  spontaneous  ignition  in  a  slab  whose  infinite  part 
is  stable. 

In  order  to  find  critical  values  of  A,  we  used  regular  continuation  in 
that  parameter  to  generate  initial  guesses  for  Newton's  method.  The  minimum 
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In  summary,  the  calculations  we  have  presented  confirm  the  usefulness  of 
the  asymptotic  boundary  conditions  described  in  [11]  and  provide  a  practical 
example  of  their  implementation. 


sendix:  Proof  of  Theorem  2.9 


We  need  to  establish  the  conditions  of  [11],  assumption  6.6.  (here,  we 
take  <5^  =  •)  Note  that  the  combinations  s(x,p;  Fu ( , 

s(x,p;  F  (u  ))(l  -  0  )  are  easily  represented  in  terms  of  the  eigenfunction 


expansion  -  name 


ly,  if  4>  =  l  c.w  , 


(A.  1  ) 


a)  s(x,p;  Fu(um))qJ>  =  \  e 

ReX^<  0 


*£(x-p) 


CiwX,'  x  >  P; 


b)  s(x,p;  FyCujXl-Qj*  =  l  e 

e 

ReXJl>0 


X^Cx-p) 


cCv  x  <  p* 


From  (2.1  Of)  we  then  have; 


(A. 2) 


-X_ (x-p) 

a)  «s(x,p;  F  (u  ))Q»  <  ae  ,  x  >  p; 

u  00  °° 


A+(x-p) 

b)  II S ( x , p ;  F  (u  )  J  (I-Q  )H  <  ae  ,  x  4  p. 

CO  "  ao  * 


..  .  ^  max..  2...  ,  max.  i,  ,  i 

We  now  estimate  llR(w  )  -  R(w  )  II  where  iw  II  <  6,  w 

x  x 

(2.4b)  and  the  definition  of  the  norm  we  have: 


ilRfw1  )  -  R(w2)l 


(A. 3) 


=  U  f  ( Uqo  +  w  ,y)  - 


fu(U»'*)W2  "  f(u«  +  W2'X}  +  fu(u-'X)w2"Wl(a)' 


We  note  that  Sobolev's  inequality  ( Agmon  [1])  implies  that: 


(A. 4)  |w*l  <  2y  HW2 II  w  {ii)  <  2Y«  . 


We  have  (suppressing  the  y  in  the  argument  of  f): 
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f(u  +  wS  -  f ( u  +  w2)  -  f  (u  )(w’  -  w2)  = 
00  2  °°2  u  00  2  2 


Jq  dt  /o  dS  fuu^U»  +  S(W2  +  t(W2  "  w2}))^w2  "  W2^w2  +  t(W2 


Taking  absolute  values  and  integrating  yields: 


Uf(u  +  w^)  -  f(u  +  w2)  -  f  (u  )(w!l  -  w2)B 
°°2  »  2  u°°22 

„  K1  1  1  2 1  ..  1  2  „ 

<  ~  2  W2  +  W2f  'iW2  ‘  Vl  (fl> 

•  dy6  «W l  -  W^I^(Q)  . 


Similarly  we  have: 

— -  (f(u  +  wS  -  f(u  +  w2)  -  f  (u  )(w^  -  w2)) 

3y.  °°  2  ®  2  u  00  2  2 

■'l 

=  (f  (u  +  wS  -  f  (u  +  w2)  -  f  ( u  )  ( w ^  -  w2)) 
y  00  2  y  oo  2  uy  00  2  2  J 

l  x  i 


+  (f  +  w,Hw!l)  -  f  (u  +  w^)  -  f  (u  )(w^  -  w't)  ) 

vu“  2  2  yi  u  “  2  yj  u  00  2  2  yi 

+  (u  )  (f  (u  +  w^)  -  f  (u  +  wj?)  -  f  (u  )(w^  -  w2) }  . 

*  V  '  u  «  2  U°°  2  uu  00  2  2  ' 

l 

The  first  and  third  terms  of  the  right-hand  side  of  the  expression  above  are 

easily  estimated  in  precisely  the  same  fashion  as  used  to  produce  (A. 5).  It 

is  only  necessary  to  replace  lfuu!  by  I f Uuv  ^  for  t*ie  ^rst  berm  and  by 

u<-  yi 

| (u  )  I  If  |  for  the  third.  For  the  second  we  have: 

1  OO '  y  I  <  UUU 
1 


1  2, 


f  (u  +  w^!)(w^)  -  f  (u  +  w2)(w2)  -  f  (u  )(w^  -  w2) 

u°“  2  2y,  u  00  2  2  y,  u"  2  2y. 

l  l  l 

=  ^0  dt  fuu^U»  +  W2  +  “  W2^(W2  "  W2^W2  +  t(W2  "  W2^y. 


+  dt  ^0  ds  fuu(u-  +  S(W2  +  t(W2  "  w2))^w2  ‘  W2^y.^W2  +  t(w2  " 
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Taking  absolute  values  and  integrating  the  expression  above  yields,  in 
combination  with  estimates  of  the  other  terms: 


l  "aTT  (f(u»  +  vl]  ~  f(u»  +  ~  fu(u«.)(w2  "  vS)li 


i  3yi 


2'  L  («) 


(A. 6) 


<  K2.2Y6«w]  -  w2.^,  +  K3.2Y6«w^  -  w2.^ 

K  K 

+  2L‘2Y6,W2  "  W2“w ,(Q)  +  I1*^6  l  MW2  -  W2)y.llL,(«) 
2  l  l  2 


From  (A.3),  (A. 5)  and  {A. 6)  we  conclude: 


(A. 7) 


ttRfw1)  -  R(w2)i  <  2y<S  (K  +  K2  +  K3)Hw1  -  w2li. 


which  implies: 


(A. 8) 


m*X#/X  dP  S(x,p)Qoo(R(w1)  -  R(w2))  -  J“  dp  S(  x  ,p)  ( I-Q^  )  (  R(  w^  ) 
<  2y6o(k1  +  k2  +  k3)(1-  +  j~)  mJx  llw1  -  w2» 


max  1  2 

<  aw  -  w  H  • 

x 


Similarly  we  have: 


(A. 9) 


BRfw1)!!  =  If (Uoo  +  wj)  -  fu(uM)w2  -  f<uj*w  (Q)  • 

2  _ 


Using  the  previous  estimates  (with  w  =  0)  we  find 


(A. 10)  llRtw1)#  <  Y62(K1  +  K2  +  K3>, 


which  implies: 

(A. 11)  m®X  «/X  dp /*  dp  S(x,p)QooR(w1)  -  /"  dp  S(  x ,  p)  ( I-Q^  )  R(  w 1  )  I!  < 


R(w2))» 
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